1. Field of the Invention
The present invention relates generally to the field of geosciences and more particularly to seismic data processing. Specifically the invention relates to a method to extract the time-lapsed changes in 3D seismic data sets collected over a production period to integrate with production data and assist in understanding and managing the extraction of oil and gas from reservoirs or the injection of other fluids into the reservoirs.
2. History of the Related Art
In the oil and gas industry, seismic surveys are carried out in order to provide subsurface images so that accumulations of hydrocarbons or other fluids might be identified. In a seismic survey, one or several sources emit elastic waves in the form of pressure or ground motion modulation from specific locations (wavefield), at or below the land or sea surface or in a borehole. This wavefield propagates away from the source(s) through the subsurface. Along with this propagation, a fraction of the incident wavefield is reflected from the fraction of the global heterogeneities in the elastic material properties of the subsurface (such as acoustic impedance). This excitation by the incident wavefield generates a reflected wavefield from the heterogeneities, which manifests as pressure, particle motion or some derived quantities and can be detected and recorded at the surface or in a borehole at a number of receiver locations.
Processing of the measurements is undertaken so as to construct a 3D image of the sub-surface. Repeated surveys at selected time intervals (days, months, years) allow observation of the changes in, over or under a given reservoir across the time interval—e.g. before oil or gas production starts and after some period of production or injection and to compare the results of measurements. This is called 4D seismic and involves comparing 3D seismic surveys carried out at different time instances. The aim is to observe changes in the state of the formations and fluids consequent upon production of hydrocarbons from or the injection of fluids into a reservoir. Proper detection of the changes and proper identification of the effects, factors and processes requires specialised acquisition techniques and data processing steps.
Such a technique applied to detect 4D changes is hereafter referred to as 3D warping. The data within the seismic data sets is first realigned or conditioned to compensate for variations in acquisition (or non-repeatability of seismic surveys) and changes in velocity in the sub-surface. The standard technique makes use of cross-correlation between different vintages in selected windows. Such a window is a time interval representing a portion of a trace. The window is set across traces for correlation and thus should contain all the 4D effects. One problem with these correlation-based approaches is the size of the correlation window. If the window used for correlation is too large, the accuracy of correlation is likely to be affected: indeed, the correlation value will then depend not only on differences between the survey at the point being considered, but also on other effects, apart from the points being considered. If the window used for correlation is too small, correlation is likely to be severely affected by noise and non-repeatability of the surveys, including changes due to the effects whose observation is desired.
In EP 1 865 340 to the Applicant, and incorporated herein by reference, the evolution of an oil reservoir in the process of producing is carried out by jointly inverting for the changes in the propagation times and seismic amplitudes of a seismic wavelet along propagation paths in the ground. Inverting allows us to back filter, in effect, deriving the original from the solution. A base survey of the reservoir is provided, with a set of seismic traces at a first time T associated to a first velocity field Vb; a monitor survey of the reservoir is provided, the monitor survey being taken at a second time T+ΔT, with a set of seismic traces associated to the same positions as in the base survey; the monitor survey is associated to a second velocity field Vm. For a set of samples i in the base survey, one computes over the samples of the set the sum S of a norm of the difference between                the amplitude bi of the seismic trace in the base survey at each sample i and        the sum of the amplitude mi′ of the seismic trace at a time-corresponding i′ in the monitor survey and the amplitude due to the reflectivity change local to the time-corresponding sample i′ induced by the difference between the first velocity field Vb and the second velocity field Vm; the time-corresponding sample i′ being shifted in time by a time-shift derived from the velocity changes along the propagation path from the surface to time-corresponding sample i′. This sum is minimised to derive the velocity changes from the base survey to the monitor survey and thus characterise the evolution of the oil reservoir.        
This analysis is based on the fact that changes in the reservoir, due to exploitation, will cause changes to the petrophysical properties of the rock and therefore to the seismic velocity field. Practically, oil will be substituted by gas or water and/or the fluid pressure will change, causing changes in saturation, porosity, permeability and pressure, and consequently changes in velocity. Changes within the reservoir may also change the stress and strain state of the surrounding rocks, further causing changes in their velocities. These changes to velocity will produce time shifts in the seismic expression of underlying reflectors and associated changes in reflectivity, causing a change in the local wavefield. By using an inversion technique, for every point in the 3D volume, an estimate of the 4D changes having occurred in the time lapse between collection of the base and monitor surveys is provided. It is therefore possible to deduce a field of 4D velocity changes without having to proceed with cross correlation of the traces.
Although the 4D inversion problem seems rather easy to formulate as the minimisation of a difference between base and monitor seismic data, it is an ill-posed problem that has multiple solutions: for instance zero-mean velocity changes map into zero time-shift. Moreover the inversion becomes even more highly non-linear for fields that induce subsidence and have potentially large time shift.
In EP 1 865 340, the crucial step is in minimising the sum. Essentially this is an optimisation problem which requires minimising of the objective function or cost function over all choices of variables i.e. velocity changes that satisfy the modelled constraints. In warping, the cost function can typically be derived as
                    C        =                              ∑                          i              =              1                        N                    ⁢                                    (                                                b                  ⁡                                      (                                          t                      i                                        )                                                  -                                  m                  ⁡                                      (                                                                  t                        i                                            -                                                                        t                          s                                                ⁢                                                                              ∑                                                          k                              =                              1                                                        i                                                    ⁢                                                                                                                    Δ                                ⁢                                                                                                                                  ⁢                                V                                                            V                                                        ⁢                                                          (                              k                              )                                                                                                                                            )                                                  +                                                      (                                                                  w                        ⁡                                                  (                          t                          )                                                                    *                                                                        Δ                          ⁢                                                                                                          ⁢                          V                                                V                                            ⁢                                              (                        t                        )                                                              )                                    i                                            )                        2                                              (        1        )            where b and m are respectively the base and the monitor trace, ts is the sampling rate of the seismic data,
      Δ    ⁢                  ⁢    V    Vis the relative velocity 4D change, w is the seismic wavelet and * denotes the convolution between the wavelet and the relative velocity change to model the 4D amplitude change.
However, as in almost any inverse problem, this cost function does not go identically to zero. In fact the forward model used for this inversion, is just an approximation of the vertical propagation that, although good, implies some assumptions and therefore a residual still exists. This is a major disadvantage of the technique proposed in EP 1 865 340, as it provides a highly unstable solution. Moreover the seismic data are affected by noise and this is seen as an additional residual after the inversion process.
A well known method to deal with approximated forward models and noise is by adding a regularisation term to the cost function. Other techniques are known, but regularisation imposes constraints on the results and thus avoids overfitting of the data and the noise. In effect we constrain the solution to the point that it does not need to be minimised. The cost function can be considered as the seismic mismatch together with the regularisation:
                    C        =                                            ∑                              i                =                1                            N                        ⁢                                          (                                                      b                    ⁡                                          (                                              t                        i                                            )                                                        -                                      m                    ⁡                                          (                                                                        t                          i                                                -                                                                              t                            s                                                    ⁢                                                                                    ∑                                                              k                                =                                1                                                            i                                                        ⁢                                                                                                                            Δ                                  ⁢                                                                                                                                          ⁢                                  V                                                                V                                                            ⁢                                                              (                                k                                )                                                                                                                                                        )                                                        +                                                            (                                                                        w                          ⁡                                                      (                            t                            )                                                                          *                                                                              Δ                            ⁢                                                                                                                  ⁢                            V                                                    V                                                ⁢                                                  (                          t                          )                                                                    )                                        i                                                  )                            2                                +                      λ            ⁢                                          ∑                                  i                  =                  1                                N                            ⁢                                                f                  ⁡                                      (                                                                                            Δ                          ⁢                                                                                                          ⁢                          V                                                V                                            ⁢                                              (                                                  t                          i                                                )                                                              )                                                  .                                                                        (        2        )            
The regularisation weight λ expresses a trade off between modelling the 4D changes from seismic and imposing constraints on the solutions. There are many forms of regularisation using any function ƒ of the relative velocity change.
In order to select the optimal regularisation, a number of steps are performed on the base and monitor seismic survey data:                1. Select a number of locations or seismic traces representative of the seismic quality;        2. Warp the data in these locations for different regularisation weights;        3. Cross-plot the cost function terms i.e. seismic mismatch against the regularisation term to provide a regularisation weight-map;        4. From a range of valid solutions determined from the plot, select the best solution to the regularisation value according to the available production history and geological information of the reservoir.        5. Interpolate the optimal regularisation values across the whole common seismic survey area to obtain the time lapsed seismic image between the base and monitor surveys.        
FIG. 1 illustrates such a cross-plot, showing four distinct regions. At A, there is no solution as the warping is trapped in a local minimum and does not converge. At B, moving to stronger regularisations, we obtain under-regularised solutions where perfect fitting of the seismic occurs but the solution is non-physical. Note this is the region of solutions for the minimisation of the sum approach proposed in EP 1 865 340. At C, we have the optimal balancing between seismic fitting and regularisation and this is the area that has to be further investigated. At D, over regularised solutions where the time shift is zero everywhere, the warped trace does not change from the monitor and the difference between warped and base is a constant.
At optimal balancing (area C), a regularisation value can be selected by using production data and geological information, that gives some hints about the expected 4D anomaly. It is therefore seen that the choice of regularisation method and regularisation weight turns out to be the most critical parameter to achieve geologically plausible solutions. Once the regularisation parameter is estimated in every selected location, then a 2D map of optimal parameters is obtained by interpolation and applied to the full data sets and a 4D inverted signal is then graphed to illustrate velocity changes and the evolution of the reservoir between the base and monitor surveys.
This process of determining the optimal regularisation weight is time intensive. While the step of warping data at the selected locations for different regularisation weights (Step 2) can be achieved on a trace by trace basis and does not require inversion of the 3D data set as a whole, the time to construct the cross plot and determine a regularisation is one to two orders of magnitude greater than the time required to run the inversion on the whole survey. Consequently, this process can take months to complete and thus the results are of historical use only since further changes in the reservoir will have taken place during this period. Such a delay limits their usefulness by a reservoir engineer who needs the information as soon as possible to take decisions on field management.